% BL_OutpuGap_REStat_2 - Figures 2 and 3, and Table 1
%
% Replication files for the paper "Measuring the output gap with large datasets"
% by Matteo Barigozzi and Matteo Luciani, forthcoming on the Review of economics and statistics.
% 

% Written by Matteo Luciani, matteoluciani@yahoo.it

clear all; close all; clc; 

ML_graph_options
stampa=0; if stampa==1; disp('Results will be printed'); end

load RestatSetting                                                          % loads benchmark settings
v2struct(Setting);                                                          % "unpack" benchmark settings
[Y, Label, Name, Dates, ~, CBO] = ML_ReadDataHaver(Filename,trans,out);     % Load data in levels
period=[40 20 8 4]; tickname={'10Y','5Y','2Y','1Y'}; wj=2*pi./period;       % Parameters for plotting spectra
ID=[1 40 43 69 75 72 73 82]; nid=length(ID);                                % Variables in Table 1
lgF={'$\tau_t$','$c_{1t}$','$c_{2t}$','$c_{3t}$','$c_{4t}$','$c_{5t}$'};

nomefile='Restat_00bw_';
LS={'k-','k--','k:','k-.'};

y=ML_diff(Y);                                                               % Data in 1st Differences
[T,N]=size(y);                                                              % size of the panel
[yy, my, sy]=ML_Standardize(y);                                             % Standardize
TT=(1:T+1)';                                                                % time trend
X=NaN(T+1,N); bt=X; b=zeros(N,2);                                           % preallocates variables for detrending
J=BL_TestLinearTrend(y); J(TR1)=1;                                          % Identify variables to be detrended
[X(:,J),bt(:,J),b(J,:)]=ML_detrend(Y(:,J));                                 % Detrend variables to be detrended
X(:,~J)=Y(:,~J)-repmat(mean(Y(:,~J)),T+1,1);                                % Demean variables not to be detrended
bt(:,~J)=repmat(mean(Y(:,~J)),T+1,1);                                       % ---------------------
b(~J,1)=mean(Y(:,~J))';b(~J,2)=0;                                           % ---------------------
if GDO==1 % --------------------------------------------------------------- % Restrictions for GDO
    b(1:2,2)=mean(b(1:2,2));                                                % same slope, different constant
    bt(:,1:2)=repmat(b(1:2,1)',T+1,1)+TT*b(1:2,2)';                         % same linear trend
    X(:,1:2)=(Y(:,1:2)-bt(:,1:2));                                          % Detrended GDP GDI
    sy(1:2)=repmat(mean(sy(1:2)),1,2);                                      % same std for standardization
end        % -------------------------------------------------------------- %
Z=X./repmat(sy,size(X,1),1);                                                % BLL - Standardization  

for qq=2:5; d=qq-1;
    [f,lambda]=ML_efactors2(yy,qq,2);                                       % estimate factor loadings, aka DFM on \Delta y_t
    if GDO==1; lambda(1:2,:)=repmat(mean(lambda(1:2,:)),2,1); end           % same loadings for GDO    
    F=Z*lambda/N;                                                           % Factors in levels as in BLL
    [A0,v,AL]=ML_VAR(F,p,1);                                                % Estimate VAR    
    xi=Z-F*lambda';                                                         % idiosyncratic component        
    Z2=(Y-repmat(b(:,1)',size(Y,1),1))./repmat(sy,size(Y,1),1);             % BL standardizations: (Y_t - a), i.e. approx. centered levels, divided by sy
    b1=b./repmat(sy',1,2); b1(:,1)=0;                                       % divide slope by sy, constant=0;      
    [A,Q,Lambda,R,F00,P0s,mu,type2]=ML_DynamicFactorSS_GDO_TV...            % State-space representation
        (AL,lambda(:,1:qq),v,qq,s,xi,F,rr,I0,I1,b1,A0,TV,star);             % --------------------------
    
    [xitT,~,~,~,~,~,~,A,Lambda,R,Q,b2]=ML_DynamicFactorEM_GDO_TV...         % EM-Algorithm
        (Z2,F00,P0s,A,Lambda,R,Q,s,qq,d,p,mu,type2,model,maxiter,tresh,cc,GDO);  % ------------    
    T2=size(Z2,1); start=1;
    BT=TT*b2';                                                                  % ML estimates of linera trend
    isTV=find(type2==10)+max(p*qq,qq*(s+1));                                      % identifies TV coefficients
    for ii=1:length(isTV) % --------------------------------------------------- % TV slopes or means
        id=find(Lambda(:,isTV(ii)));
        BT(:,id)=xitT(:,isTV(ii))*ones(1,length(id));
    end                   % --------------------------------------------------- %
    FF=xitT(:,1:qq); FF1=[]; for ss=0:s;FF1=cat(2,FF1,FF(s-ss+1:end-ss,:));end   % Store ML Factors
    lambda2=NaN(N,qq,s+1);                                                  % Store ML loadings
    for ss=1:s+1; lambda2(:,:,ss)=Lambda(:,(ss-1)*qq+1:ss*qq); end          % -----------------
    L=reshape(lambda2,N,qq*(s+1));                                          % -----------------
    
    T3=length(FF1(cc:end,:));
    [V, D]= eig(cov(FF(cc:end,:)));                                             % eigenvalue and eigenvenctors of VCV
    [~, t2]=sort(diag(D),'descend'); V=V(:,t2); D=D(t2,t2);                     % sort eigenvalues and eigenvectors
    FFt = FF(cc:end,:)*V(:,1:qq-d);                                              % TREND Factorss
    FFc = FF(cc:end,:)*V(:,qq-d+1:qq);                                            % CYCLE Factors
    
    %%% =================== %%%
    %%%  Common Components  %%%
    %%% =================== %%%
    SY=repmat(sy,T3,1);
    MY=repmat(b(:,1)',T3,1);
    start2=start+s+cc-1;
    Y2=Y(start2:end,:);
    chi=(FF1(cc:end,:)*L'+BT(start2:end,:)).*SY+MY;                         % common component
    chint=FF1(cc:end,:)*L'.*SY+MY;                                          % common component without linear trend
    zeta=Y2-chi;                                                            % idiosyncratic component
        
    TableXi(qq-1,:)=var(zeta(:,ID));   % ---------------------------------- % TABLE 1 
    
    Dates2=Dates(start2:end); Dates2d=Dates2(2:end); Dates4q=Dates2(5:end); % New Dates for plots
    j0=find(year(Dates2)==1970,1);                                          % Strarting point for all graphs    
    colore={[0 0.45 0.74],[0.64 0.08 0.18],[0.85 0.33 0.1],...              % colors for gaps blue, red, orange
        [0.93 0.69 0.13],[0 0 0],[0.47 0.67 0.19],[0.49 0.18 0.56]};        % yellow,black, green, purple    
    qsp=[num2str(qq) num2str(s) num2str(p)];
    DD=Dates2(j0:end);
    
    if qq==4
        ZZ0=[FFt(j0+1:end,:) FFc(j0+1:end,[3 1 2])];
        [s1, omega]=ML_Spectrum(ML_diff([FFt FFc(:,[3 1 2])]));
    else
        [s1, omega]=ML_Spectrum(ML_diff([FFt FFc(:,end:-1:1)]));
    end
    
    %%% ====================== %%%
    %%%  Figure 2 upper plots  %%%
    %%% ====================== %%%
    if qq<5
        axes('Parent',figure,'FontSize',12); ML_FigureSize, hold on;
        for jj=1:size(s1,2); pl1(jj)=plot(omega,s1(:,jj),LS{jj},'linewidth',2); end
        hold off; axis tight; box on;
        xlim=get(gca,'xlim');  ylim=get(gca,'ylim');
        axis([xlim 0 ylim(2)+ylim(2)/20])
        set(gca,'Xtick', wj,'xticklabel',tickname,'yticklabel',num2str(get(gca,'ytick')','%.1f'));
        %     gridxy(get(gca,'xtick'),get(gca,'ytick'),'color',[.8 .8 .8],'linewidth',1)
        lgd =legend(pl1,lgF(1:qq),'interpreter','latex'); lgd.FontSize = 14;
        if stampa==1; print('-depsc','-painters','-r600',[nomefile 'SpettroFF_' num2str(qq) '.eps']); end
    end
end

%%% ========= %%% 
%%%  Table 1  %%%
%%% ========= %%%
disp(' ')
disp('Idiosyncratic variances')
disp(TableXi)
disp(' ')


J=[1 3 4 7];
ZZ0(:,4)=-ZZ0(:,4);
ZZ=ML_Standardize(ZZ0);
WW=ML_Standardize(ML_diff(ZZ0));
SizeXY=ML_SizeXY(DD,ZZ,1,2); SizeXY2=ML_SizeXY(DD(2:end),WW,1,2); % -----------------

%%% ========== %%%
%%%  Figure 3  %%%
%%% ========== %%%
for ii=1:4
    pl=ML_TimeSeriesUSbw(ZZ(:,ii),DD,[],SizeXY,LS(1));
    lgd = legend(pl,lgF(ii),'location','Best','interpreter','latex'); lgd.FontSize = 14;
    if stampa==1; print('-depsc','-painters','-r600',[nomefile 'FF' num2str(ii) '.eps']); end    
    pl=ML_TimeSeriesUSbw(WW(:,ii),DD(2:end),[],SizeXY2,LS(1));
    lgd = legend(pl,lgF(ii),'location','Best','interpreter','latex'); lgd.FontSize = 14;
    if stampa==1; print('-depsc','-painters','-r600',[nomefile 'dFF' num2str(ii) '.eps']); end
end     
    

ID2=[1 75 40 43 41];
[s1, omega]=ML_Spectrum(ML_Standardize(y(:,ID2)));
xlim=ML_minmax(omega);
ylim=[0 ML_min(s1,2)*1.025];

for ii=1:2    
    axes('Parent',figure,'FontSize',12); ML_FigureSize, hold on;
    pl1=plot(omega,s1(:,ii),'k','linewidth',1.5);
    hold off; axis tight; box on;    
    axis([xlim 0 ylim(2)+ylim(2)/20])
    set(gca,'Xtick', wj,'xticklabel',tickname,'yticklabel',num2str(get(gca,'ytick')','%.1f'));
%     gridxy(get(gca,'xtick'),get(gca,'ytick'),'color',[.8 .8 .8],'linewidth',1)    
    if stampa==1; print('-depsc','-painters','-r600',[nomefile 'SpettroY_' num2str(ID2(ii)) '.eps']); end
    title(Name{ID2(ii)},'fontsize',14,'fontweight','bold')
end


axes('Parent',figure,'FontSize',12); ML_FigureSize, hold on;
pl1=plot(omega,s1(:,3),'k','linewidth',1.5);
pl2=plot(omega,s1(:,4),'k--','linewidth',1.5);
pl3=plot(omega,s1(:,5),'k-.','linewidth',1.5);
hold off; axis tight; box on;
axis([xlim 0 ylim(2)+ylim(2)/20])
set(gca,'Xtick', wj,'xticklabel',tickname,'yticklabel',num2str(get(gca,'ytick')','%.1f'));
% gridxy(get(gca,'xtick'),get(gca,'ytick'),'color',[.8 .8 .8],'linewidth',1)
legend([pl1 pl2 pl3],{'Total','excluding Food & Energy','Energy'},'fontsize',12)
if stampa==1; print('-depsc','-painters','-r600',[nomefile 'SpettroY_' num2str(ID2(3)) '.eps']); end




